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Abstract 

We determine mass and mixing angles of r/ and if states using Nf = 2 + 1 + 1 
Wilson twisted mass lattice QCD. We describe how those flavour singlet states 
need to be treated in this lattice formulation. Results are presented for three values 
of the lattice spacing, a = 0.061 fm, a = 0.078 fm and a = 0.086 fm, with light 
quark masses corresponding to values of the charged pion mass in a range of 230 
to 500 MeV and fixed bare strange and charm quark mass values. We obtain = 
557(15) (45) MeV (first error statistical, second systematic) and (j) = 44(5)° for a 
single mixing angle in the quark flavour basis, = —10(5)° in the octet-singlet basis. 



1 Introduction 

From experiment it is known that the masses of the nine light pseudo-scalar 
mesons show an interesting pattern. Taking the quark model point of view, the 
three lightest mesons, the pions, contain only the two lightest quark flavours, 
the up- and down-quarks. The pion triplet has a mass of M n m 140 MeV. For 
the other six, the strange quark contributes also, and hence they are heavier. 
In contrast to what one might expect, five of them, the four kaons and the rj 
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meson, have roughly equal mass around 500 to 600 MeV, while the last one, 
the rf meson, is much heavier, with a mass of about 1 GeV. On the QCD 
level, the reason for this pattern is thought to be the breaking of the Ua(1) 
symmetry by quantum effects. The rj' meson is, even in a world with three 
massless quarks, not a Goldstone boson. 

In this paper, we present the first lattice study of rj and rj' meson masses using 
mass degenerate up, down and heavier, non-degenerate strange and charm 
dynamical quark flavours. The lattice QCD formulation is the Wilson twisted 
mass formulation [1] with Nf = 2 + 1 + 1 dynamical quark flavours. This will 
not only allow a study of the dependence of the 77, 77' masses on the light quark 
mass value, but also an investigation of the charm quark contribution to both 
of these states. Moreover, the r\ c meson can be studied in the unitary case in 
principle. 

rj and rj states are difficult to treat in lattice QCD, because fermionic discon- 
nected contributions appear and cannot be ignored. This is why the amount 
of available results for these states from lattice QCD is rather limited in range 
of pion mass values as well as values of the lattice spacing. For recent lat- 
tice studies in Nf = 2 + 1 flavour QCD see [2,3,4,5]. Some of these we will 
discuss in more detail later and compare with our results. A previous lattice 
study with Wilson twisted mass fermions has been performed with Nf = 2 
dynamical quark flavours [6]. In this work we have developed particular noise 
reduction techniques, which we also used to obtain the results presented here 
for the N f = 2 + 1 + 1 case. 

The paper is organised as follows. In section 2 we introduce the lattice QCD 
framework we are using, the Wilson twisted mass formalism, followed by a 
discussion of how we deal with flavour singlet pseudo-scalar mesons in this 
framework in section 3. In section 4 we present our results and discuss them 
in the following final section. More details on our analysis procedure and data 
tables can be found in the appendix. 



2 Lattice Action 

The lattice QCD action for Nf = 2 + 1 + 1 Wilson twisted mass fermions reads 



For the gauge action S g we use the Iwasaki gauge action [8]. The twisted mass 
Dirac operator for the light - i.e. up/down quark - doublet reads [1] 



S = S g + xeDeXe + Xh D h Xh ■ 



(1) 



D £ = D w + m + ifia 5 r 3 



(2) 



2 



ensemble 


P 






ans 


L/a 


A^conf 


N s 


N b 


A30.32 


1.90 


0.0030 


0.150 


0.190 


32 


1367 


24 


5 


A40.24 


1.90 


0.0040 


0.150 


0.190 


24 


2630 


32 


10 


A40.32 


1.90 


0.0040 


0.150 


0.190 


32 


863 


24 


4 


A60.24 


1.90 


0.0060 


0.150 


0.190 


24 


1251 


32 


5 


A80.24 


1.90 


0.0080 


0.150 


0.190 


24 


2449 


32 


10 


A100.24 


1.90 


0.0100 


0.150 


0.190 


24 


2493 


32 


10 


A80.24s 


1.90 


0.0080 


0.150 


0.197 


24 


2517 


32 


10 


A100.24s 


1.90 


0.0100 


0.150 


0.197 


24 


2312 


32 


10 


B25.32 


1.95 


0.0025 


0.135 


0.170 


32 


1484 


24 


5 


B35.32 


1.95 


0.0035 


0.135 


0.170 


32 


1251 


24 


5 


B55.32 


1.95 


0.0055 


0.135 


0.170 


32 


1545 


24 


5 


B75.32 


1.95 


0.0075 


0.135 


0.170 


32 


922 


24 


4 


B85.24 


1.95 


0.0085 


0.135 


0.170 


24 


573 


32 


2 


D15.48 


2.10 


0.0015 


0.120 


0.1385 


48 


1045 


24 


10 


D45.32sc 


2.10 


0.0045 


0.0937 


0.1077 


32 


1887 


24 


10 



Table 1 

The ensembles used in this investigation. The notation of ref. [7] is used for labeling 
the ensembles. In addition we give the number of configurations iV con f, the number 
of stochastic samples N s for all ensembles and the bootstrap block length iV&. 

and for the strange/charm doublet [9] 

D h = D w + m + i^a^r 1 + fi S r 3 , (3) 

where Dy/ is the Wilson Dirac operator. The value of m was tuned to its 
critical value m crit as discussed in refs. [10,7] in order to realise automatic 
O(a) improvement at maximal twist [11]. This way of 0(a) improvement was 
first shown to work in practice in refs. [12,13] in the quenched approximation, 
and later for Nf = 2 in ref. [14] (see also ref. [15]). For a review see ref. [16]. 
Also for Nf = 2 + 1 + 1 we have strong indications that the scaling properties 
of Wilson twisted mass lattice QCD at maximal twist are favourable [7]. For 
a discussion on how to determine kaon and D-meson masses see ref. [17]. The 
bare twisted masses \x a and fj,s are related to the bare strange and charm quark 
masses via the relation 

m CyS = fi a ± Zfi s (4) 

where Z = (Zp/Zs) denotes the ratio of pseudo-scalar and scalar renormali- 
sation constants Zp and Z$. Quark fields in the twisted basis are denoted by 
Xt,h and in the physical basis by ipe,h- They are related via the axial rotations 

^ = e^ 3 l\ t , ^ = xt e iw75r3/4 , 

With automatic 0(a) improvement being the biggest advantage of twisted 
mass lattice QCD (tmQCD) at maximal twist, the downside is that flavour 
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p 


1.90 


1.95 


2.10 




5.231(38) 


5.710(41) 


7.538(58) 



Table 2 

Values of rj/a for the three 0- values. 

(and parity) symmetry is broken at finite values of the lattice spacing. This 
was theoretically and numerically shown to affect mainly the mass value of the 
neutral pion [14,18,19], however, in the case of Nf = 2+1 + 1 dynamical quarks, 
it implies the complication of mixing between strange and charm quarks. 

We use gauge configurations as produced by the European Twisted Mass Col- 
laboration (ETMC) with action described above [7,17,20]. The details of the 
configurations are described in ref. [7] and the ensembles used in this investi- 
gation are summarised in table 1: we use ensembles denoted with A, B and 
D with values of the lattice spacing a A = 0.0863(4)fm, a B = 0.0779(4)fm and 
a D = 0.0607(2)fm, corresponding to /3 A = 1.90, (3 B = 1.95 and (3 D = 2.10, 
respectively [20]. The physical volumes are with only a few exceptions larger 
than 3 fm. In the table we also compile the number of investigated gauge 
configurations and the number of stochastic samples per gauge configuration 
used to estimate the disconnected contributions. 

Throughout this paper we will use the Sommer parameter tq [21] to investigate 
the scaling of our results. We use the values of rj/a extrapolated to the chiral 
limit at each f3- value separately. For f3 — 1.90 and /3 — 1.95 we use the values 
quoted in ref. [7]. For (5 = 2.10 we did the extrapolation by ourselves. All values 
for rj/a are shown in table 2. For setting the physical scale we could use the 
results of ref. [20], where a chiral fit to data for /pg and mps and the physical 
value of f w was used to set the scale. However, the fit in ref. [20] includes only 
two data points at j3 — 2.10, and hence, is rather preliminary. Therefore, we 
prefer to use a value of r = 0.45(2) fm in this paper to set the scale. The 
5% error covers the statistical uncertainty and spread quoted in ref. [20] and 
allows room for systematic uncertainties. As soon as an update of the scale 
setting becomes available the results in this paper can be updated accordingly. 
For fixing the light and strange quark masses to their physical values we will 
use the experimental values of M n o = 135 MeV and M K o = 498 MeV. We use 
the masses of the neutral mesons to reduce uncertainties from the fact that 
we do not include electromagnetic effects in our simulation. See also ref. [22] 
and references therein for a discussion. 

For every /3-value the bare values of a/i a and afi$ were kept fixed, and only 
the value of the light twisted mass parameter a/j,£ is varied. The kaon masses 
measured on these ensembles are close to the physical value with a deviation 
of up to 10% [7] in particular for the A and D ensembles. The D-meson mass 
values have a large uncertainty, but they are also close to physical [7]. We will 
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discuss this point later in more detail. 

For the ensembles A80.24 and A100.24 we have additional ensembles with re- 
tuned values of a/j, a and a^s denoted by A80.24s and A100.24s (see table 1), 
which reproduce the physical kaon mass value more accurately than the origi- 
nal A80.24 and A100.24 ensembles (see figure 4). We will use these ensembles 
to estimate the strange quark mass dependence of the 77 mass. 



3 Flavour singlet pseudo-scalar mesons 



In order to compute masses of pseudo-scalar flavour singlet mesons we have 
to include light, strange and charm contributions to build the appropriate 
correlation functions. In the light sector, one appropriate operator is given 
by [6] 

-^(ipuhb^u + ^Pdhb^d) ->■ -j={-XuXu + XdXd) = O e , (6) 

in the physical and the twisted basis, respectively. With twisted mass fermions 
we have to work with doublets of quarks, hence, in the strange and charm 
sector the corresponding operators read 

The sign in 1 ± r 3 in the physical basis distinguishes the charm and strange 
quark contribution. As a consequence, working in the twisted basis we need 
to compute correlation functions of the following interpolating operators 

O c = Z(xJl5Xc-Xsil5Xs)/2- (XsXc + XcXs)/2, , . 

O s = Z(xsH5Xs ~ XdlhXc)/^ ~ (XsXc + XcXs)/ 1 } ■ 

Note that the sum of pseudo-scalar and scalar contributions appears with 
the ratio of renormalisation factors Z, which needs to be taken into ac- 
count properly. The renormalistaion factors for eq. (7) are the non-singlet 
Zp and Z$. While singlet and non-singlet Zp are identical, the singlet and 
non-singlet Z$ differ at two loop order [23] in perturbation theory. Defining 
O = (Oi, O c , O s ) , the correlation function matrix is given by 

C(t)=^(O(x)^O(0)) . (9) 

X 

However, masses are independent of Z as well as the choice of basis, so, for the 
purpose of determining masses we can proceed as follows: starting from the 
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bilinears in eqs. (6) and (8) we change the operator basis via an appropriate 
rotation matrix 1Z 



c n (t) = (no(t) <g> no(o)) = nc(t)n T , 



(10) 



where 

^VOiOt Vo e S h Vo e P h \ 

C n = Vs h o e Vs h s h Vs h v h , (11) 
VtPhOt Vr h s h WhVhJ 

is again a symmetric, real and positive definite correlation matrix with 

V h = {xdlbXc ~ X^75Xs)/2 , S h = Z~ l (xsXc + XcXs)/2 (12) 

and 7]xy denoting the corresponding correlation function. The rotation matrix 
1Z is given by 

A \ 



K 



n I L 

+73 -73/ 



(13) 



Now we can drop the factor Z' 1 , which appears only as a constant scaling 
factor for Sh- We denote the corresponding correlation matrix by Cr. This 
disentanglement of scalar and pseudo-scalar contributions greatly reduces the 
number of terms required for each element of the correlator matrix. 



Solving the generalised eigenvalue problem [24,25,26] 

Cn(t) V (n) (t,t ) = A< B >(Mo) Cn(to) v {n) (t,t ) , 



(14) 



and taking into account the periodic boundary conditions for a meson, we can 
determine the effective masses by solving 



m (")(T-t) 



A (w) (Mo) _ 
A< n )(t + 1, t ) e - m(n) (* +1 ) + e - m(n) (r-(t+i)) 



(15) 



for , where n counts the eigenvalues. The state with the lowest mass should 
correspond to the r\ and the second state to the r( meson. Alternatively, we 
use a factorising fit of the form 



A A , 

2m( n ) 



exp(-m (n) t) + exp(-m w (T - *)) 



(16) 



to the correlation matrix matrix C. For this we either first rotate Cn back 
taking the factor Z into account, see below, or we directly construct C taking 
Z into account. The amplitudes A q>n correspond to (0\qq\n) with n = r),r)', ... 
and q — £,s, c. 
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3.1 Flavour Content and Mixing 



One might first of all be interested in the quark flavour content of a given 
state in order to compare to phenomenology. From the components 77o?i,2 of 
the eigenvectors rf n ^ defined above, we can reconstruct the flavour contents 

(n) 
C £,s,c 

(13), we reconstruct cj l ] c from 



of the states. Since we have changed the basis according to eqs. (11) and 



(n) _ _J_ f „(«K 

c^-j^i-Z-^ + V^)/^ (17) 

4 B) = ^o(-^ B) -^ B) )/v^ 

with normalisation 

A/» = v / ri n) ) 2 + (^- 1 ^ n) ) 2 + (^ n) ) 2 . 

At this point the ratio Z is again required. The flavour non-singlet renormal- 
isation factors have been evaluated non-perturbatively for our situation [27]. 
Another way to get access to these ratios of renormalisation constants is avail- 
able: we require that unphysical amplitudes (such as the connected correlator 
from source sTs to sink cTc in the physical basis), which are formally order a 
contributions in the twisted mass formulation, be minimised. 

For instance, for 535 this procedure gives Z = 0.70, while for D15 we ob- 
tain Z = 0.75. For comparison the preliminary RI-MOM flavour non-singlet 
value [28] for f3 = 1.95 is 0.700(8) and at = 2.1 is 0.737(14). This close 
agreement indicates that we have a reliable estimate of Z. Moreover, it turns 
out that the consequence of an error in this ratio on our mixing angles is 
minimal, see the discussion in section 4. 

The mixing between r] and rf is usually expressed in terms of mixing angles 
in an appropriate basis. Here we use the quark flavour basis considering only 
light and strange quarks. With 

\Ve) = -^(\uu) + \dd)) , \r} s ) = \ss) 

one arrives at 

^ /cos0 -sin<A _ f\r)e)\ 
ysincf) cos0 J \\r) s )J 
For a detailed discussion including the charm quark see for instance ref . [29] . 

On the lattice, however, we have to work with the amplitudes A q>n defined 
above. Following refs. [2,5], we first rotate the matrix C-n (t) back to the original 
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form C (t) in eq. (9) in the way prescribed above (including Z). The amplitudes 
A q n in eq. (16) of a factorising fit to the rotated correlation matrix are then 
directly related to the mixing angles via [29] 



'A e:V A S A = ff e cos fa -/ s sin0 s \ 
A e w A s „> fs cos <f) s j ' 



where we ignored the charm contribution. Hence, the mixing angles <pe and <p s 
can be extracted from 

tan0 £ = ^-, tan0 s = -^^-, (19) 

where the renormalisation constants cancel in the ratio. Following the RBC / 
UKQCD and the Hadron Spectrum collaborations we also define a common 
angle (representing the geometric mean of <pg and <f> s ) 

tan 2 (0) = -^^, (20) 

inspired by arguments that <pe and <p s should actually agree [29,30], see also 
ref. [5]. We do not include the charm quark in the discussion here, because it 
turns out that its contribution to rj and rj' is negligible. 



4 Results 



We have computed all contractions needed for building the correlation ma- 
trix of eq. (11). For the connected contributions, we used stochastic time-slice 
sources (the so called "one-end-trick" [31]). For the disconnected contribu- 
tions, we used stochastic volume sources with complex Gaussian noise [31]. 
As discussed in ref. [6] one can estimate the light disconnected contributions 
very efficiently using the identity 

D' 1 ~ D d 1 = -2i^Dt 75 D- 1 . 

For the heavy sector such a simple relation does not exist, but we can use 
the so called hopping parameter variance reduction, which relies on the same 
equality as in the mass degenerate two flavour case (see ref. [31] and references 
therein) 

D h l = B- BHB + B(HBf - B{HBf + D h 1 (HB) 4 

with Dh = (1 + HB)A, B = 1/A and H the two flavour hopping matrix. We 
use 24 stochastic volume sources per gauge configuration in both the heavy 
and the light sector. 
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ensemble 


oMpg 


aM K 


aM v 


aM v > 


A30.32 


0.12374(27) 


0.25150(29) 


0.286(15) 


0.49(6) 


A40.24 


0.14517(39) 


0.25884(43) 


0.281(18) 


0.39(6) 


A40.32 


0.14174(26) 


0.25666(23) 


0.281(11) 


0.49(9) 


A60.24 


0.17340(39) 


0.26695(52) 


0.290(7) 


0.59(9) 


A80.24 


0.19888(37) 


0.27706(61) 


0.302(8) 


0.60(8) 


A100.24 


0.22097(40) 


0.28807(34) 


0.315(11) 


0.50(7) 


A80.24s 


0.19870(50) 


0.25503(33) 


0.270(9) 


0.54(9) 


A100.24s 


0.22149(39) 


0.26490(74) 


0.280(5) 


0.66(12) 


B25.32 


0.10768(30) 


0.21240(50) 


0.234(10) 


0.50(8) 


B35.32 


0.12445(29) 


0.21840(28) 


0.237(9) 


0.59(9) 


B55.32 


0.15503(48) 


0.22799(34) 


0.249(14) 


0.60(10) 


B75.32 


0.18121(24) 


0.23753(32) 


0.253(13) 


0.44(7) 


B85.24 


0.19373(58) 


0.24392(59) 


0.260(12) 


0.51(12) 


D15.48 


0.07004(19) 


0.16897(85) 


0.201(9) 


0.38(8) 


D45.32sc 


0.07981(30) 


0.17570(84) 


0.192(15) 


0.30(4) 



Table 3 

Results of aM n , aM n i for all ensembles and the corresponding values for the charged 
pion mass Mpg and the kaon mass Mr. Most of the kaon and pion mass values have 
been published already in ref [7], however, we have recomputed the pion mass values 
with our statistics. 



ensemble 






L 


4>f 




4> F 


A30.32 


57(15) 


49(13) 


56(15) 


36(9) 


53(13) 


46(11) 


A40.24 


39(16) 


41(11) 


39(16) 


43(11) 


40(12) 


41(13) 


A40.32 


56(18) 


33(11) 


54(18) 


33(11) 


44(11) 


44(11) 


A60.24 


66(15) 


29(13) 


66(15) 


30(13) 


49(8) 


49(8) 


A80.24 


73(7) 


20(10) 


73(7) 


24(10) 


48(6) 


50(6) 


A100.24 


62(10) 


37(10) 


62(10) 


39(9) 


50(5) 


51(5) 


A80.24s 


59(16) 


36(16) 


59(16) 


41(16) 


48(12) 


50(12) 


A100.24s 


75(13) 


20(16) 


76(12) 


25(16) 


50(7) 


54(7) 


B25.32 


74(13) 


25(10) 


73(13) 


28(10) 


52(6) 


53(6) 


B35.32 


67(14) 


19(7) 


63(15) 


22(8) 


42(6) 


41(7) 


B55.32 


74(11) 


14(6) 


73(11) 


18(7) 


43(6) 


46(6) 


B75.32 


50(14) 


52(15) 


50(14) 


52(15) 


51(11) 


51(11) 


B85.24 


60(16) 


47(20) 


61(16) 


46(19) 


54(13) 


54(13) 


D15.48 


59(14) 


25(9) 


57(14) 


29(9) 


41(10) 


42(10) 


D45.32sc 


60(19) 


57(15) 


59(19) 


58(15) 


59(15) 


58(15) 



Table 4 

Results for the mixing angles in eqs. (19) and (20) from a 4 x 4-correlation function 
matrix using local (L) and fuzzed (F) operators for all ensembles 

Also, we use both local and fuzzed operators to enlarge our correlation matrix 
by a factor two. In addition to the interpolating operator quoted in eqs. (6) 
and (7), one could also consider the 7-matrix combination 27075, which would 
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5 10 15 5 10 15 

t/a t/a 

(a) (b) 

Fig. 1. (a) Effective masses in lattice units determined from solving the GEVP for 
a 3 x 3 matrix with to /a = 1 for ensemble B25.32. (b) the same, but for a 6 x 6 
matrix. For comparison also the fit results (see text) for the two lowest states are 
shown (cf. table 3). 

increase the correlation matrix by another factor of two. However, the corre- 
sponding correlation functions turn out to be too noisy to give any further 
improvement at this stage. The number of gauge configurations iV con f investi- 
gated per ensemble and the number of stochastic sources per gauge configu- 
ration N s are summarised in table 1. 



The details of our GEVP and fitting procedures are explained in appendix A, 
together with fit ranges and \ 2 values. We point out here that our experience 
is that it is rather important to sample as many gauge configurations as pos- 
sible. Autocorrelation plays a significant role for the rf state (in agreement to 
what we found in the Nf = 2 case [6]), whereas the r\ is less affected. Errors are 
always computed using a bootstrap procedure with 1000 bootstrap samples 
and blocking. Block sizes have been chosen such that the distance in between 
blocks in units of HMC trajectories of length r = 1 (see ref. [32] for a de- 
scription of the HMC algorithm used) is at least 20. Therefore, the statistical 
significance of the number of gauge configurations N con{ in table 1 depends on 
the ensemble, i.e. the number of trajectories in between two consecutive gauge 
configurations. We could not work with longer blocks, because the number of 
blocks becomes too small. The number of configurations per block Nf, is also 
given in table 1. 



4-1 Extraction of Masses 



In the left panel of figure 1 we show the effective masses determined from solv- 
ing the GEVP for ensemble 525.32 from a 3 x 3 matrix with local operators 
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only. For the purpose of this plot we kept to/ a — 1 fixed. One observes that the 
ground state is very well determined and it can be extracted from a plateau 
fit. The second state, i.e. the rj', is much more noisy and a mass determination 
is questionable, at least from a 3 x 3 matrix. Enlarging the matrix size sig- 
nificantly reduces the contributions of excited states to the lowest states and, 
due to smaller statistical errors at smaller t values, a determination becomes 
possible, which one can see for a 6 x 6 (including fuzzing) matrix in the right 
panel of figure 1. Our final results and errors for the masses are indicated in 
figure 1 by the horizontal bands. We do not determine them from a constant 
fit to the plateau, but by a three state cosh fit to the eigenvalues, possibly 
with larger to values, see appendix A. As shown, the procedure gives very 
good agreement with a plateau fit for the 77, but slightly lower values for the 
7/. This indicates no n- negligible systematic uncertainties in the extraction of 
the 7/ mass value. 

The third state appears to be in the region where one would expect the r] c mass 
value, however, the signal is lost at t/a — 5, which makes a reliable determi- 
nation unfeasible. Note that the two plots in figure 1 are rather independent 
of the particular ensemble chosen. 

To gain further confidence in our identification of the 77 and 7/ states, we 
also determine the flavour content of the two states as explained above. As 
an example we show in the left panel of figure 2 the flavour content of the 
77 for ensemble 525.32. It becomes evident that - as one would expect from 
phenomenology - the 77 has a dominant strange quark content, while the 7/', 
shown in the right panel of figure 2 is dominated by light quarks. Note that 
we do not include a gluonic operator in our analysis, so we discuss only the 
relative quark content. For both, the charm contribution is compatible with 
zero. 

In figure 3 we show the masses of the 77 (filled symbols) and 7/ (open symbols) 
mesons for the various ensembles we used as a function of the squared pion 
mass, everything in units of r®. All the masses have been determined from 
solving the GEVP for a 6 x 6 matrix, the details are described in appendix A. 
The results have been independently cross-checked using a factorising fit, and 
corroborated. We have collected all the values for aM n and aM v > together with 
kaon and pion mass values in table 3 with statistical errors only. Note that 
the values for czMk and aMpg are published for most of the ensembles already 
in ref. [7]. But we have recomputed the pion mass values with our statistics, 
and added ensembles as compared to ref. [7] . The methods for computing the 
kaon mass in Nf = 2 + 1 + 1 Wilson twisted mass lattice QCD are described 
in ref. [17]. 

It is clear from the figure that the 77 meson mass can be extracted with high 
precision, while the 7/' meson mass is more noisy. The former can be under- 
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Fig. 2. (a) Squared flavour content of rj for B25.32 from 3 x 3-matrix using local 
operators only, (b) Squared flavour content of rf for B25.32. 
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Fig. 3. (a) r\ (filled symbols) and rf (open symbols) masses in units of chirally 
extrapolated tq (listed in table 2) as a function of (ro-Mps) 2 . 

stood, as in the SU(3) symmetric limit the rj meson is a flavour octet with all 
disconnected contributions vanishing, while the rj' is the flavour singlet with 
non-vanishing disconnected contribution^]. 

The results displayed in figure 3 have been obtained using the bare values 



We thank Martin Savage for a useful discussion on this point. 
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Fig. 4. (a) The kaon mass in units of vq and (b) r^M v as a function of (roMps) 2 
also for all A ensembles. The dotted and solid curves in (a) represent the fitted qk 
in eq. (23) and the shifted gx (see text), respectively. 

of an a and afis as used for the production of the ensembles. Those values, 
however, did not lead to the physical values of, e.g., the kaon and D- meson 
masses [7,17]. We show the kaon mass as a function of the squared pion mass 
in the left panel of figure 4 for all A ensembles. It is clearly evident that the re- 
tuned ensembles (A80.24s and A100.24s) have values for the kaon mass closer 
to the physical one (see also e.g. ref. [7,20]). In the right panel of figure 4 we 
show tqM v also for all A ensembles (including v480.24s and y4100.24s) and 
the same pattern as for the kaon mass is observed. Furthermore, the physical 
strange and charm quark mass values differ among the A, B and D ensembles. 
Hence, figure 3 is not yet conclusive with regards to the size of lattice artifacts 
and the extrapolation to the physical point. What we can recognise is that 
the light quark mass dependence in the r\ appears to be rather weak. 

While the results for the rj mass in figure 3 show a consistent picture over all 
lattice spacings and light quark mass values - keeping the differences in the 
strange and charm quark masses in mind - the rj' mass shows large fluctuations 
and no consistent picture. We attribute this to two observations: firstly, due 
to the large noise in the rj' state, the extracted masses are for most of the 
ensembles only an upper bound for M^, because a plateau is hardly visible. 
Therefore, we think the values of the rf mass have very likely a non-negligible 
systematic uncertainty stemming from the fact that the signal for the rf is 
lost in noise at rather small t-values. We are currently investigating a solver 
for the GEVP using a singular value decomposition, which is better suited to 
deal with the noise. First results are encouraging. 

Secondly, the strange (and charm) quark mass values are different among the 
three values of the lattice spacing. The strange quark mass value has a strong 
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Fig. 5. (a) roMrj as a function of (a/ro) 2 for the ensembles ^460.24, B55.32 and 
D45.32sc. (b) Values of r^M^ as a function of (roMps) 2 as explained in the text. The 
continuum extrapolated values r^M^ in (a) and the experimental value roM V!exp , 
obtained with tq = 0.45(2) fm in (b) are horizontally displaced for legibility. 

influence on the rf meson mass, as we learn from A80.24, A80s, A100.24 and 
Al00.24s, see figure 3. 



4-2 Scaling Artifacts and Strange Quark Mass Dependence of M n 



For M v the statistical uncertainty is sufficiently small to allow for a meaningful 
scaling test. For this we need to compare M v at the three different values of 
the lattice spacing for fixed values of for instance r M K , r M D , r M PS and the 
physical volume. From volume and the charm quark mass value we expect only 
little influence given our uncertainties and hence, we are going to disregard 
effects from slightly different physical volumes at the different /3-values and 
the differences in the charm quark mass in the following. 

As we do not have simulations at the three values of the lattice spacing with 
matched values of r M^, we have to perform an interpolation in Mk- For 
this procedure we have to rely on two pairs of ensembles, namely (A80.24, 
A8Q.2As) and (A100.24, A100.24s). The two ensembles within a pair differ in 
the values of a/z CT and a/i<s, whereas afii is identical. We can use these ensembles 
to estimate the derivative of M 2 with respect to M| and use this estimate 
to correct for the mismatch in r M^. By using this estimate of the derivative 
for all ensembles - not only A80.24 and ^4100.24 - we neglect the dependence 
of Djj on the lattice spacing and the light and charm quark masses. 

In more detail, we treat the masses of the //-meson and the kaon like in chiral 
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perturbation theory as functions M 2 = M 2 [Mp S , M£] and define the dimen- 
sionless quantity 



d(aM v ) 2 
d(aM K ) 2 



(21) 



Next we make the approximation that D v is independent of the quark mass 
values fj,£, /i a , us an d (3. Its value is actually equal within errors when estimated 
from A80.24 and ,480.24s or from A100.24 and ,4100.24s and on average we 
obtain D v = 1.60(18). 

Now we use this value of D v to correct the three ensembles ,460, B55 and D45 
- which have approximately equal values of r M PS ~ 0.9 - to a common value 
of r M K ~ 1.34 using 

(r M ri ) 2 = (r M v ) 2 + D v • A K , 

where Ak is the difference in the squared kaon mass values (in units of ro). We 
plot the resulting tqM^ values for the three ensembles A60, B55 and DA5 as a 
function of (a/r ) 2 in the left panel of figure 5. The data are compatible with 
a constant continuum extrapolation roM^° st = 1.480(34), which we indicate 
by the horizontal line. We can also attempt a linear extrapolation, which is 
also shown in the figure, leading to r M"J? n ° = 1.61(14). The difference in 
between the two extrapolated values 

r AM^° = 0.13(13) (22) 

gives us an estimate on the systematic uncertainty to be expected from the 
continuum extrapolation. We cannot repeat this analysis for more values of 
r Mps, as we have currently only two D ensembles analysed. We will therefore 
quote a 8% relative error from AM^ /^! f° r our mass estimates, which 
we believe is a conservative figure. 



Nevertheless, in order to obtain a more complete picture, we now attempt to 
correct all our ensembles for the slightly miss-tuned value of Mk- For this we 
adopt the following procedure: first we perform a linear fit 

9K[(r M PS ) 2 , a,b}=a + b- (r M PS ) 2 (23) 

to the values of (r M K ) 2 for the A-ensembles (without A80.2As and A100.24s). 
We obtain a = 1.492(6) and b = 0.571(8). This curve is then shifted to a value 
a = 1.238(6) such that (r M£ xp ) 2 = g K [(r M n ) 2 ,a,b\. We denote this new 
function by g K = g K [(r Mp S ) 2 ,a,b}. The two curves corresponding to Qk and 
gK are shown in the left panel of figure 4. As it turns out, the kaon masses 
for the two tuned ensembles A80.24s and Al00.24s are already very close to 
c/k- Now we can correct the remaining A ensembles as well as the B and D 
ensembles to the same line of g by computing the difference of the squared 
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kaon mass values to g 



6 K [(r M PS ) 2 ] = (r M K ) 2 [(r Mp S ) 2 ] - ~g K [(r M PS ) 2 } , 

which we can use to correct the measured M 2 [(r M PS ) 2 ] corresponding to 

(r M,) 2 [(r M PS ) 2 ] = (r M,) 2 [(r Mp S ) 2 ] + D v ■ 5 K [(r M PS ) 2 ] . (24) 

The result of this procedure is shown in the right panel of figure 5: we show 
values of r^M^ for all our ensembles as a function of (r M PS ) 2 . It is evident 
that all the data fall on a single curve within statistical uncertainties. Figure 5 
confirms that M v is not affected by large cut-off effects. Note again that we 
ignored the /ig, \i a , fj,s and j3 dependence of D v and that we cannot fully 
estimate the systematics stemming from this approximation. 

4-3 Extrapolation to the Physical Point 

As gKroM^) 2 ] = (r M^ xp ) 2 - and the strange quark mass was fixed to its 
physical value using M£ xp - we can next attempt a linear fit to all cor- 
rected (by the procedure discussed above) data points for (roM v ) 2 [(r M P o,) 2 ]. 
Using r = 0.45(2) as discussed in section 2, the fit yields r^M^ [r^M 2 ] = 
1.252(58) st at(100) sys and in physical units 

M v (M n ) = 549(33) stat (44) sys MeV , 

where the experimental mass- value of the neutral pion M n o = 135 MeV has 
been used for M w . In the SU(2) chiral limit we obtain r M° = 1.230(65) sta t(98) s 
or 

M° = 539(35) stat (43) sys MeV . 

For estimating the systematic error we used the ratio AM^/Mji as dis- 
cussed above. Note that the error on tq = 0.45(2) significantly contributes to 
the statistical errors for the results in physical units. 

As the procedure presented above relies on the assumption that D v is inde- 
pendent on fj, a , us and /3, it would be desirable to have a cross-check. 
Inspecting figure 4 once more, one observes that M n and Mk appear to have 
a similar strange quark mass dependence. This motivates to search for an 
appropriate ratio of quantities in which the strange quark mass dependence 
cancels approximately. One option is to study the ratio M V /M K , which is 
shown in the left panel of figure 6 as a function of (r M P s) 2 for all available 
ensembles. Strikingly, within errors all data points fall on the same curve, and 
in particular the points for ,480.24s and ,480.24s and A100.24 and ,4100.24s 
agree within errors, respectively. This confirms in particular that most of the 
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Fig. 6. (a) M^/Mk as a function of (ro-Mpg) 2 for all available ensembles, (b) GMO 
ratio (see text) as a function of (ro-Mpg) 2 for all available ensembles. Experimental 
values are horizontally displaced for legibility. 

strange quark mass dependence cancels in the ratio. Extrapolating all the data 
for (M v /Mk) 2 linear in (roMpg) 2 to the physical pion mass point we obtain 



{MJMk] 



1.121(26) 



which is in good agreement with the experimental value (M v /M^ 
Using the experimental value of M K o = 498 MeV we obtain 



cxp 



(25) 



1.100. 



M v = 558(13) stat (45) sys MeV, 

where the first error is statistical and the second systematical as estimated 
from the scaling violations discussed above. The latter represents a rather 
conservative estimate since some of the scale dependence might cancel in the 
ratio M^/Mk- Note that in this procedure the scale r = 0.45(2) fm enters only 
for determining the physical pion mass point. As the slope of the extrapolation 
is rather small, the statistical uncertainty in M v is smaller than for the direct 
extrapolation of (r M v ) 2 . 

Alternatively, motivated from chiral perturbation theory, we also study the 
GMO relation 



3M: = AMI 



Mi 



^ — - • (26) 
It is valid in the SU(3) symmetric case, but violated only by a few percent 
with physical values of the corresponding meson masses. Therefore, the strange 
quark mass dependence of 3M 2 /(4M|- — M 2 ) should be weak. We show the 
dimensionless ratio 3M 2 /(4M£ — M 2 ) as a function of (r Mpg) 2 in the right 
panel of figure 6. The first interesting remark on this figure is that again the 
data points for the ratio from A80.24 and ,480.24s (A100.24 and Al00.24s) - 
which differ in the bare strange quark mass - agree within errors, confirming 
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Fig. 7. Comparison of our results for M„ (filled symbols) (corrected for the mismatch 
in Mk as discussed in the text) and Mw (open symbols) in physical units for all 
three values of the lattice spacing to results from the literature (RBC/UKQCD [2], 
HSC [4], UKQCD [5]. The scale for our points was set using tq = 0.45(2) fm. 

that a large part of the strange quark mass dependence cancels in the ratio 
(like for M^/M-k). Again, all the data points fall onto one single line within 
errors, independent of the value of the lattice spacing, the strange and the 
charm quark mass. If we fit a linear function in (r Mps) 2 to the data we 
obtain 

(3M*/(4M K - Ml))^ = 0.966(48) (27) 

at the physical pion mass, which is in agreement with experiment, (3M 2 / (4M K — 
M^)) cxp = 0.925. Using the experimental values of M n o and Mro we now ob- 
tain 

M v = 559(14) stat (45) sys MeV , 

where the first error is statistical and the second systematical estimated from 
the scaling violations discussed above. 

We remark that our value for D v = 1.60(18) is in rather good agreement to 
the value 4/3 one would obtain naively from eq. (26). 

In figure 7 we show a compilation of our results for r\ masses (corrected for 
the mismatch in Mk, cf. figure 5) and rj' masses and those available in the 
literature for Nj = 2 + 1 flavour lattice QCD: in Ref. [3] r] and rj' meson masses 
have been computed using Nj = 2 + 1 flavours of overlap quarks at one value 
of the lattice spacing and large values of the pion mass. In this reference not 
enough details are given to be included in our comparison figure 7. In Ref. [2] 
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these masses have been determined using Nf = 2 + 1 flavours of domain wall 
fermions, again at a single value of the lattice spacing a ~ 0.1 fm and three 
values of the pion mass ranging from about 400 MeV to about 700 MeV. The 
corresponding data points are labeled RBC/UKQCD in figure 7. In Ref. [4] the 
Hadron Spectrum Collaboration (HSC) - also at a single value of the lattice 
spacing - and one value of the pion mass used Wilson fermions. Finally, in 
Ref. [5] data is presented for two values of the pion mass, each at a different 
value of the lattice spacing with staggered fermions, labeled as UKQCD in 
our plot. As discussed previously, we used a value of r = 0.45(2) fm to set 
the scale for our data. 

Firstly, from figure 7 it is clear that the results presented in this paper sig- 
nificantly increase the available data for rj and rf masses from lattice QCD. 
In particular, our results extend to significantly lower values of the pion mass 
than available before. And, within the statistical uncertainties, all the results 
from the different lattice formulations agree, despite the fact that the system- 
atic uncertainties are not taken into account, and that we used Nf = 2 + 1 + 1 
dynamical quark flavours. 

Secondly, figure 7 allows to compare the statistical uncertainties quoted by 
the different collaborations. As stated previously, our experience is that our 
statistical error is dominated by gauge noise, and not by stochastic noise. This 
is in particular so for the r/, where also autocorrelation times are large. For 
this reason we had to sample more than 1000 gauge configurations (cf. table 1) 
for most of the investigated ensembles, and still our statistical error might be 
slightly underestimated, as we could not make the blocks in our bootstrap 
analysis arbitrarily long. The RBC/UKQCD and HS collaboration results are 
obtained on a significantly smaller set of gauge configurations, and hence it is 
unclear to us how those errors compare to ours. 



4-4 V and rj' Mixing 

As mentioned previously, the determination of the rj and rj' mixing angles 
requires the knowledge of the ratio of renormalisation constants Z as discussed 
above. 

We determine the amplitudes A^, A^, A SjV and A SjV > from a factorising fit to 
the correlation matrix C(t) rotated to its original form. For the factorising fit 
we consider only light and strange degrees of freedom (both, local and fuzzed), 
as the charm does not contribute within errors. We then determine <pe and <p s 
from eq. (19) and <j) from eq. (20). The errors are again determined using a 
bootstrap procedure with 1000 bootstrap samples. All results from a 4 x 4 
matrix fit are collected in tab. 4. Results from a 6 x 6 matrix can be found in 
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Fig. 8. (a) mixing angle (f> L from local amplitudes of a 4 x 4-matrix (no charm quark), 
(b) mixing angle cj> F from fuzzed amplitudes of 6 x 6-matrix (including charm quark) 



the appendix. 



The results for are shown in figure 8. In the left panel we show </> in degrees 
determined from local amplitudes (L), and in the right panel from fuzzed 
amplitudes (F), both as a function of (r Mps) 2 . Both agree nicely and we do 
not observe a dependence on the lattice spacing and (roMx) 2 , at least within 
the relatively large errors. The dependence on (r M PS ) 2 is weak, the data is 
compatible with a linear extrapolation to the physical pion mass point, which 
leads to a value of 

<\> = 44(5)° , 

using a combined fit to local and fuzzed data. In the octet basis this would 
correspond to a rjg, r}\ mixing angle 9 = <p — 54.7° = —10(5)°. For this extrap- 
olation we ignored a possible strange and charm quark mass dependence of 
0. Our result for is in agreement with the results from RBC/UKQCD [2], 
HSC [4] and an old UKQCD work [33], which all quote something in between 
40° and 50°. In the recent UKQCD staggered investigation [5] a value of 34(3)° 
is favoured. Comparing to experimental and phenomenological results [29,34,5] 
we find excellent agreement to results from radiative decays and glueball mix- 
ing (~ 42°). Results from photon fusion and charm-77 production favour a 
value of ~ 33°. Note that one has to keep in mind that our determination of 
the mixing angle is likely to be affected by systematic uncertainties. 

We have checked that this determination is not affected by our uncertainty 
on the ratio of renormalisation constants Z. For this we varied Z in a large 
range from 0.4 to 1.0, and we found that this does not affect the extraction of 
the mixing angles. As an example we show the dependence of <fi L for ensemble 
D15.48 in figure 9 on Z. The variation introduced by Z is by far smaller 
than the statistical uncertainty. We observe the same for the other ensembles 
and conclude that our evaluation of the angles is not affected by systematic 
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Fig. 9. Z-dependence of the double-ratio mixing angle <j) L denned in eq. (20) for 
D15.48 using local amplitudes only. 

uncertainties stemming from the Z ratio within our statistical uncertainties. 

While figure 8 indicates that the angle behaves regularly and the linear 
extrapolation to the physical pion mass appears to be reasonable, the results 
for <pe and (fi s are diverse, see table 4. We attribute this to large statistical and 
systematic uncertainties in these quantities. We hope to be able to investigate 
4>i and 4> s further once we have improved our determination of the r{ state. 



5 Summary and Discussion 

We presented the first computation of r\ and rf meson masses from lattice QCD 
with degenerate up/down, heavier strange and even heavier charm dynamical 
quarks in the Wilson twisted mass formulation. Our results based on ETMC 
gauge configurations cover three values of the lattice spacing from ~ 0.09 fm 
down to ~ 0.06 fm and a range of pion masses from 230 MeV to 500 MeV. 
The results presented here are therefore the first covering three values of the 
lattice spacing and a range of pion masses down to 230 MeV. Kaon and D- 
meson masses are close to their physical values. We observe in both, rj and rj' 
masses, only a mild dependence on the light quark mass. The r] meson mass 
can be determined with high statistical accuracy, while the rj suffers from 
noise and large autocorrelation times. 

In our simulations the bare strange and charm quark masses were kept fixed 
for each value of the lattice spacing separately, apart from two additional A- 
ensembles. As the kaon mass values deviate from the experimental values by 
up to 10%, we applied three methods to correct for this mismatch: firstly, 
by estimating the strange quark mass dependence from the aforementioned 



so 
70 
60 
50 
40 
30 
20 
10 



21 



additional A-ensembles, secondly by considering the ratio M v /Mk and thirdly, 
by using the GMO relation eq. 26. All three methods seem to indicate that 
lattice artifacts are only weakly affecting the rj mass. As our best estimate we 
extract 

M v = 557(f5) stat (45) sys MeV 

by a weighted average over the three methods. The first error is statistical 
and the second systematic from the continuum extrapolation. The result is in 
good agreement with the experimental value for M v . 

We also determine the 77 and 7/ mixing angle <fi using amplitudes from a fac- 
torising fit model. Again, we observe little light quark mass and lattice spacing 
dependence within errors. But this time also the strange quark mass depen- 
dence is smaller than our (large) statistical uncertainty. Our best estimate for 
the mixing angle is 

= 44(5)° 

with statistical error only. In an octet basis this corresponds to 6 = —10(5)°. 
The agreement of our estimate for with most experimental and phenomeno- 
logical [29,34] and most other lattice determinations [2,4,33] is excellent. The 
value is also very close to the quadratic GMO estimate of 44.7°. But results 
from photon fusion and charm-?? production, as well as the lattice determi- 
nation from ref. [5] find a value of ~ 34°. Estimating the systematics for the 
angle becomes difficult, because the statistical uncertainties for each ensemble 
are large. 

It is evident that the rj' mass is affected by large statistical and systematical 
uncertainties. It also shows a large autocorrelation time. Hence, we currently 
cannot make a reliable estimate of the physical rj' mass value. We therefore 
plan in the future to apply additional noise reduction techniques to get a bet- 
ter signal also for the 7/ state. There are actually two promising approaches, 
the point-to-point method described in ref. [6] and an extension of the noise 
reduction trick used here only for estimating the light disconnected contribu- 
tions. Also, investigating the SVD solver further might be a promising way 
forward. 

Finally, determining the flavour singlet decay constants is of large phenomeno- 
logical interest. We are currently investigating our signals for (0\A^\rj) and 
(OlAJr/). 
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Table A.l 

Parameters of the GEVP applied to the 6x6-matrix from local and fuzzed operators. 
A Analysis Details 

In this appendix we summarise the technical details of our numerical methods 
for extracting masses and angles. 

A.l Generalised Eigenvalue Problem (GEVP) 

For the determination of masses and identification of flavour contents of the 
states we employ the variational approach in eq. (14). We extract the eigen- 
values A to) from the GEVP and determine then the masses from a fit 
to A. The errors of A are determined using a bootstrap procedure with 1000 
bootstrap samples, after the data is blocked in order to account for autocor- 
relation. An overview of the remaining GEVP and fitting parameters is listed 
in table A.l. Due to the rather large difference in the signal-to- noise ratio 
between the two lowest lying states we use two different approaches for r] and 
T)': 

(1) For the ground-state (77) we fit a single cosh in t in a region \t\,t^[ to 
our data for \ v (t,t ). The lower bound of the fit-range t\ is chosen by 
visual inspection of the effective mass plot to lie at the beginning of the 
plateau and also such that further increasing t\ does not change the value 
of the resulting mass within errors. The latter is also true for the choice 
of the starting value t[j in eq. (14). In principal choosing t larger leads 
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Table A. 2 

Values for Z from matching mixed action (Osterwalder-Seiler) to unitary approach. 

to smaller masses, though due to noise we are only able to moderately 
increase t before the error gets too large. 
(2) For the first excited state (rj 1 ) we perform a three state cosh-fit to the data 
of \ r ''(t,t ), starting from the lowest t\ possible, i.e. t\ = t ' +1 in order 
to use as many points as possible. This is necessary because for many 
ensembles there is no clear plateau in the effective masses reached before 
the signal is lost in noise. Therefore, this procedure is a major source 
of systematic error for the determination of rj' masses, at least for those 
cases where only few points are available. Only for the D-Ensembles it 
turns out to be possible to choose t$ > 1, for all other ensembles we had 
to use t[J = 1. 

The upper bound of the fit range t% is independently determined for every 
state n by the last eigenvalue (£2, to) distinguishable from noise. Note that 
the value of £2 is not very important for the fit, as eigenvalues at large t have 
typically large errors and therefore do not contribute much to the fit. 

A. 2 Factorising fit model 

For the extraction of amplitudes which are required to calculate mixing angles 
we employ the factorising fit model as detailed in eq. (16) to the matrix in 
the original twisted basis in eq. (9). We limit ourselves to n = 2 in eq. (16), 
i.e. a two state fit. As discussed before the ratio of renormalisation constant 
Z = Zp/Z s does not affect masses and angles, but only amplitudes. For 
the results of the angles we used the values summarised in table A. 2. They 
have been obtained by matching a mixed action to the unitary action and 
extrapolating to the chiral limit, see ref. [38]. These values agree well with the 
RI-MOM determination of ETMC [28,27], and the method we discussed in 
the text. 

As we have already determined the masses from the GEVP, we use those 
together with their errors (see table 3) as priors to our factorising fit. This 
stabilises the fit procedure, but we have checked that it does not affect the 
result. We always use an uncorrelated x 2 function, because the correlation 
matrix is too noisy. Table A. 3 lists the input parameters and the resulting 
uncorrelated x 2 /dof for the fits to the 4x4 (light and strange degrees of 
freedom, local and fuzzed) and 6x6 (light, strange charm, and local and 
fuzzed) correlator matrix. The results quoted in the main text for the angles 
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Table A. 3 

Parameters for the factorising fit. 

were determined using the 4x4 matrix, because the charm does not contribute. 
However, with the 6x6 matrix we obtain almost identical results. 

The value of the lower bound of the fit range t\ is chosen to be constant for 
every value of the lattice spacing, whereas the choice of the upper bound t 2 is 
determined by the requirement that x 2 /dof < 1, and that it is close to the t% 
used for GEVP. 
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